///////////////////////////////////////////////////////////////////////////////
//
// labFisStefBolt.cc - acquire data for SNS Lab Fis resonance experiment
// 
// 1st release: A. Lusiani <a.lusiani@sns.it> April 2003
//
// this program includes parts from poll.c, Part of Comedilib,
// Copyright (c) 1999,2000 David A. Schleef <ds@schleef.org>
// poll.c may be freely modified, distributed, and combined with
// other software, as long as proper attribution is given in the
// source code.
//

///////////////////////////////////////////////////////////////////////////////
//
// Modifications:
//
//

///////////////////////////////////////////////////////////////////////////////
//
// headers
//

extern "C" {
#include <stdio.h>
#include <stdlib.h>
#include <termios.h>
#include <comedilib.h>
#include <fcntl.h>
#include <unistd.h>
#include <errno.h>
#include <ctype.h>
#include <string.h>
#include <math.h>
#include <fcntl.h>
#include <time.h>
#include <sys/select.h>
}

#ifndef __CINT__
#include "TFile.h"
#include "TApplication.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TNtuple.h"
#include "TRandom.h"
#endif

#include "labFisDaq.h"
#include "aluSleep.hh"

///////////////////////////////////////////////////////////////////////////////
//
// function declarations
//

///////////////////////////////////////////////////////////////////////////////
//
// variable  declarations
//

const char* 		COMEDI_DEV = "/dev/comedi0";
int			AI_SUBDEV;
int			AO_SUBDEV;
//--- comedi sub-devices for Digital Output
int			DO_SUBDEV;

//--- buffer size for streaming data out of the driver
const unsigned int	COMEDI_BUFSIZ =	(4*1024*1024);

//--- Analog Input definitions
const unsigned int      TEMP1_AI_CHAN=1;
const unsigned int      TEMP2_AI_CHAN=2;
const unsigned int 	VBULB_AI_CHAN=5;
const unsigned int 	VSHUNT_AI_CHAN=6;

//
// input channels for asynchronous DAQ
// - A5 bulb voltage
// - A6 bulb + shunt resistor voltage
//
#define AISCAN_BEG	(5)
#define AISCAN_NUM	(2)
unsigned int aiScanList[AISCAN_NUM];

//--- Analog Input voltage range
const float 		AI_VMIN = 0;
const float 		AI_VMAX = 10;

//--- select analog input voltage reference from GROUND
const unsigned int	AI_AREF = AREF_GROUND;

//--- rate at which the inputs are sampled in Hz
unsigned int 		AISCAN_RATE = 1000;
//--- number os samples to get for each voltage
#define 		AISCAN_SAMPLES (1000)

//--- Analog Out definition
const unsigned int 	OVBULB_AO_CHAN = 0;

//--- minimum Analog Output value for smallest frequency
double 			OVBULB_MIN = 0;
//--- maximum Analog Output value for largest frequency
double			OVBULB_MAX = 8;
//--- number of separate voltages applied
const unsigned int      OVBULB_NUM = 8;

//--- select Analog Output voltage reference from GROUND
const unsigned int	AO_AREF = AREF_GROUND;

///////////////////////////////////////////////////////////////////////////////
//
// code
//

int main(int argc, char *argv[])
{
  int i;
  int rc;
  comedi_t *cifh;
  comedi_cmd cmd;
  lsampl_t maxDataAi;
  lsampl_t maxDataAo;
  double minVoltsAi, maxVoltsAi;
  double minVoltsAo, maxVoltsAo;
  unsigned int rangetAi;
  comedi_range* rangepAi;
  unsigned int rangetAo;
  comedi_range* rangepAo;
  //--- buffer for ni_pcimio Analog Input driver (16-bit unsigned samples)
  sampl_t buf[AISCAN_SAMPLES*AISCAN_NUM];
  
  //--- endpoints of acquisition are valid numbers
  rc = comedi_set_global_oor_behavior(COMEDI_OOR_NUMBER);

  cifh = comedi_open(COMEDI_DEV);
  if(!cifh){
    perror(COMEDI_DEV);
    exit(1);
  }

  //--- set comedi device streaming buffer size
  rc = comedi_set_buffer_size(cifh, AI_SUBDEV, COMEDI_BUFSIZ);
  if (rc<0) {
    fprintf(stderr, "comedi_set_buffer_size: %s\n",strerror(errno));
    exit(1);
  } else {
    rc = comedi_get_buffer_size(cifh, AI_SUBDEV);
    if (rc<0) {
      perror("comedi_get_buffer_size");
    } else {
      fprintf(stderr, "comedi buffer size for AI set to 0x%x bytes\n", rc);
    }
  }

  //--- find which subdevice offers Analog Input
  rc = comedi_find_subdevice_by_type(cifh, COMEDI_SUBD_AI, 0);
  if (rc<0){
    comedi_perror("comedi_find_device AI");
    exit(1);
  }
  AI_SUBDEV = rc;

  //--- find which subdevice offers Analog Output
  rc = comedi_find_subdevice_by_type(cifh, COMEDI_SUBD_AO, 0);
  if (rc<0){
    comedi_perror("comedi_find_device AO");
    exit(1);
  }
  AO_SUBDEV = rc;

  //--- find AI range for -10 Volt .. 10 Volt
  rc = comedi_find_range(cifh, AI_SUBDEV,
			 VBULB_AI_CHAN, UNIT_volt, AI_VMIN, AI_VMAX);
  if (rc<0){
    comedi_perror("comedi_find_range");
    exit(1);
  }
  rangetAi = rc;

  //--- get ptr to range struct
  rangepAi = comedi_get_range(cifh, AI_SUBDEV, VBULB_AI_CHAN, rangetAi);
  //--- get maximum possible sample for AI
  maxDataAi = comedi_get_maxdata(cifh, AI_SUBDEV, VBULB_AI_CHAN);
  minVoltsAi = comedi_to_phys(0, rangepAi, maxDataAi);
  maxVoltsAi = comedi_to_phys(maxDataAi, rangepAi, maxDataAi);

  fprintf(stderr, "Comedi offers AI range n. %d, %f .. %f\n",
	 rangetAi, minVoltsAi, maxVoltsAi);

  //--- build channel list for DAQ scan of microphone
  for(i=0;i<AISCAN_NUM;i++){
    aiScanList[i] = CR_PACK(AISCAN_BEG+i, rangetAi, AI_AREF);
  }

  //--- find AO range for -10 Volt .. 10 Volt
  rc = comedi_find_range(cifh, AO_SUBDEV,
			 OVBULB_AO_CHAN, UNIT_volt,
			 OVBULB_MIN, OVBULB_MAX);
  if (rc<0){
    comedi_perror("comedi_find_range");
    exit(1);
  }
  rangetAo = rc;

  //--- get ptr to range struct
  rangepAo = comedi_get_range(cifh, AO_SUBDEV, OVBULB_AO_CHAN, rangetAo);
  //--- get maximum possible sample for AO
  maxDataAo = comedi_get_maxdata(cifh, AO_SUBDEV, OVBULB_AO_CHAN);
  minVoltsAo = comedi_to_phys(0, rangepAo, maxDataAo);
  maxVoltsAo = comedi_to_phys(maxDataAo, rangepAo, maxDataAo);

  fprintf(stderr, "Comedi offers AO range n. %d, %f .. %f\n",
	 rangetAo, minVoltsAo, maxVoltsAo);

  //--- prepare AI scan for microphone's input
  labFisScanPrepare(cifh, AI_SUBDEV, &cmd,
		    TRIG_NOW, AISCAN_RATE, TRIG_COUNT, AISCAN_SAMPLES,
		    aiScanList, AISCAN_NUM);
  
  //
  // initialize Root
  //
  // code inspired from:
  // - test/hsimple.cxx of Root 3.05.3 package
  //

  //
  // Create a new ROOT binary machine independent file.
  // This file is now becoming the current directory.
  //
  TFile hfile("labFisStefBolt.root", "RECREATE",
	      "Lab Fis: Legge di Stefan-Boltzmann");

  //--- allocate Root Tree
  TTree *tree = new TTree("nt", "Legge di Stefan-Boltzmann");

  //--- Root ntuple definition: see Root/tutorials/cernbuild.C

  //
  // ntuple variables
  // - set voltage across bulb
  // - measured voltage across bulb
  // - voltage across shunt resistor
  // - temperatures from probes: T1, T2 
  //

  // - time in 1/100 s
  // Int_t hsecs;
  // tree->Branch("hsecs", &hsecs, "hsecs/I");

  Float_t ntVbSet;
  Float_t ntVb;
  Float_t ntVs;
  Float_t ntT1;
  Float_t ntT2;
  //--- add branches, each with a single leaf
  tree->Branch("vbSet", &ntVbSet, "vbSet/F");
  tree->Branch("vb", &ntVb, "vb/F");
  tree->Branch("vs", &ntVs, "vs/F");
  tree->Branch("t1", &ntT1, "t1/F");
  tree->Branch("t2", &ntT2, "t2/F");

#define HS_TOTSAMPLES ((OVBULB_NUM+1)*AISCAN_SAMPLES)
#define HS_TOTCHANNS  (HS_TOTSAMPLES/20)
#define HS_XMIN  (-0.5/AISCAN_RATE)
#define HS_XMAX  ((static_cast<double>(HS_TOTSAMPLES)-0.5)/AISCAN_RATE)
#define HS_YMIN  (-10)
#define HS_YMAX  (10)

  //--- voltage set across bulb
  TProfile *hsVbulbSet =
    new TProfile("hsVbulbSet", "Voltage set across bulb",
		 HS_TOTCHANNS, HS_XMIN, HS_XMAX, HS_YMIN, HS_YMAX);
  hsVbulbSet->SetXTitle("seconds");
  hsVbulbSet->SetYTitle("Volts");

  //--- voltage measured across bulb
  TProfile *hsVbulb =
    new TProfile("hsVbulb", "Voltage measured across bulb",
		 HS_TOTCHANNS, HS_XMIN, HS_XMAX, HS_YMIN, HS_YMAX);
  hsVbulb->SetXTitle("seconds");
  hsVbulb->SetYTitle("Volts");

  //--- voltage measured across shunt resistor
  TProfile *hsVshunt =
    new TProfile("hsVshunt", "Voltage across shunt resistor",
		 HS_TOTCHANNS, HS_XMIN, HS_XMAX, HS_YMIN, HS_YMAX);
  hsVshunt->SetXTitle("seconds");
  hsVshunt->SetYTitle("Volts");

  //--- ask for non-blocking I/O on comedi
  fcntl(comedi_fileno(cifh), F_SETFL, O_NONBLOCK);

  //--- allocate object to setup terminal input
  labFisTermIoHandler tioh(stdin);
  //--- setup terminal for immediate keypress input with no echo
  tioh.SetImmedNoecho();

  {
    //--- set bulb voltage at 0 Volts
    lsampl_t data = comedi_from_phys(0, rangepAo, maxDataAo);
    fprintf(stderr, "Set bulb voltage at 0 Volts (0x%04x, %4d)\n",
	    data, data-0x8000);
    rc = comedi_data_write(cifh, AO_SUBDEV, OVBULB_AO_CHAN,
			   rangetAo, AO_AREF, data);
    if (rc<0){
      comedi_perror("comedi_write_data");
      exit(1);
    }
  }

  //--- wait some seconds
  aluSleep(4.);
  // system("echo -n \007");
  
  //
  // read temperatures
  //
  const int TEMPLOOPNUM = 20;
  double t1Sum = 0, t2Sum = 0;
  for(int i = 0; i < TEMPLOOPNUM; i++) {
    double tVal;
    
    fprintf(stderr,"Temperatures:");
    for(unsigned int iChan=TEMP1_AI_CHAN; iChan<=TEMP2_AI_CHAN; iChan++) {
      lsampl_t data;
      
      rc = comedi_data_read(cifh, AI_SUBDEV, iChan, rangetAi, AI_AREF, &data);
      if (rc<0){
	comedi_perror("comedi_read_data");
	exit(1);
      }
      tVal = comedi_to_phys(data, rangepAi, maxDataAi);
      fprintf(stderr," %8.4f", tVal);
      switch (iChan) {
      case TEMP1_AI_CHAN: t1Sum += tVal; break;
      case TEMP2_AI_CHAN: t2Sum += tVal; break;
      }
    }
    fprintf(stderr,"\n");
  }
  ntT1 = t1Sum / TEMPLOOPNUM;
  ntT2 = t2Sum / TEMPLOOPNUM;

  //
  // set bulb voltage to several voltages
  // for each voltage, acquire voltage and current for some time
  //
  for (unsigned int iVbulb = 0; iVbulb <= OVBULB_NUM; iVbulb++) {
    //--- set ntuple variable to requested bulb voltage
    ntVbSet = OVBULB_MIN+
      static_cast<double>(iVbulb)/OVBULB_NUM*(OVBULB_MAX-OVBULB_MIN);
    
    {
      //--- set bulb voltage
      lsampl_t data = comedi_from_phys(ntVbSet, rangepAo, maxDataAo);
      fprintf(stderr,
	      "Setting bulb voltage, step %3d, %2.2fV (0x%04x, %4d)\n",
	      iVbulb, ntVbSet, data, data-0x8000);
      rc = comedi_data_write(cifh, AO_SUBDEV, OVBULB_AO_CHAN,
			     rangetAo, AO_AREF, data);
      if (rc<0){
	comedi_perror("comedi_write_data");
	exit(1);
      }
    }

    // fprintf(stderr, "start\n");
    //--- wait some seconds to get to equilibrium
    aluSleep(2.);
    // fprintf(stderr, "stop\n");

    //--- start comedi Analog Input scan
    rc = comedi_command(cifh, &cmd);
    fprintf(stderr,"Start DAQ, rc= %d\n", rc);
    if (rc<0) {
      fprintf(stderr,"errno=%d\n",errno);
      comedi_perror("comedi_command");
      exit(1);
    }

    //--- init sample number
    unsigned int scansTot = 0;
    
    //--- DAQ loop
    for(;;) {
      //--- read available data
      rc = labFisScanRead(cifh, AI_SUBDEV, &buf[0], sizeof(buf),
			  static_cast<long int>(1e6), 10);
      
      //--- end endless loop if end of data reached
      if (rc == 0) {
	fprintf(stderr, "EOF, ending scan after reading %6d samples\n",
		scansTot);
	break;
      }
      
      //--- read available data
      if (rc > 0) {
	unsigned int i;
	
	unsigned int len = rc/sizeof(sampl_t);
	for(i=0; i < len; i++, scansTot++) {
	  unsigned int scanNum = scansTot / AISCAN_NUM;
	  double vtime = 
	    static_cast<double>(iVbulb)*
	    AISCAN_SAMPLES/AISCAN_RATE;
	  vtime += static_cast<double>(scanNum)/AISCAN_RATE;
	  switch(scansTot % AISCAN_NUM) {
	  case 0:
	    //--- modulo #0 samples: bulb voltage
	    //--- fill histogram of read voltages vs. time in seconds
	    ntVb = comedi_to_phys(buf[i], rangepAi, maxDataAi);
	    hsVbulb->Fill(vtime, ntVb);
	    break;
	  case 1:
	    //--- modulo #1 samples: bulb + shunt voltage
	    //--- fill histogram of read voltages vs. time in seconds
	    ntVs = comedi_to_phys(buf[i], rangepAi, maxDataAi);
	    ntVs -= ntVb;
	    hsVshunt->Fill(vtime, ntVs);

	    //--- store this event ntuple after setting last DAQ variable
	    tree->Fill();
	    hsVbulbSet->Fill(vtime, ntVbSet);
	    break;
	  }

	}
	continue;
      }
      
      if (-rc == LFDAQ_STDINRDY) {
	unsigned char ch = fgetc(stdin);
	if (ch == 25) {
	  fprintf(stderr,"  --- <control-Y> received, quitting ... --- \n");
	  break;
	} else {
	  fprintf(stderr,"  --- Press <control-Y> to quit --- \n");
	}
	continue;
      }
      
      if (-rc == LFDAQ_TIMEOUT) {
	fprintf(stderr,"Timeout waiting for data \n");
	break;
      }
      
      break;
    }
  }

  {
    //--- leave bulb voltage at 0 Volts
    Float_t ntVbSetLast = 0;
    lsampl_t data = comedi_from_phys(ntVbSetLast, rangepAo, maxDataAo);
    fprintf(stderr,
	    "\nLeaving bulb voltage at %2.2fV (0x%04x, %4d)\n",
	    static_cast<double>(ntVbSetLast), data, data-0x8000);
    rc = comedi_data_write(cifh, AO_SUBDEV, OVBULB_AO_CHAN,
			   rangetAo, AO_AREF, data);
    if (rc<0){
      comedi_perror("comedi_write_data");
      exit(1);
    }
  }

  //--- Save all objects in this file
  hfile.Write();

  //--- Close the file.
  hfile.Close();

  tioh.Reset();

  fprintf(stderr,
          "Please select \"File->Quit Root\" on Root window to quit\n");

  //--- open saved file for reading
  TFile hfileR("labFisStefBolt.root", "READ");

  //--- recover histograms
  TProfile *hsVbulbSetR = (TProfile *)hfileR.Get("hsVbulbSet");
  hsVbulbSetR->SetMaximum(8);
  hsVbulbSetR->SetMinimum(-0.5);
  TProfile *hsVbulbR = (TProfile *)hfileR.Get("hsVbulb");
  hsVbulbR->SetMaximum(8);
  hsVbulbR->SetMinimum(-0.5);
  TProfile *hsVshuntR = (TProfile *)hfileR.Get("hsVshunt");
  hsVshuntR->SetMaximum(1.5);
  hsVshuntR->SetMinimum(-0.1);

  //--- init an object to get Root GUI on a canvas
  TApplication theApp("App", &argc, argv);

  //--- init an object to get a graphical canvas
  TCanvas *canvas = new TCanvas("canvas", "Bulb voltage vs. time", 700, 700);
  canvas->Divide(1,2);

  //--- draw histograms
  canvas->cd(1);
  hsVbulbR->Draw("hist C");
  canvas->cd(2);
  hsVshuntR->Draw("hist C");

  //--- get the canvas updated
  canvas->Update();

  //--- recover pointer to Root tree (ntuple) "nt"
  // TTree *treeR = (TTree *)hfileR.Get("nt");

  //--- start root GUI: "Exit ROOT" will return to the next statement
  theApp.Run(kTRUE);

  //--- Close the file.
  hfile.Close();

  return 0;
}
